function [TQ_tab] = mt_inv_table_depth(TQ_tab,stationlist,...
    vertchannel_station_list_amp,vertchannel_station_list_phase,...
    eastchannel_station_list_amp, eastchannel_station_list_phase, ...
    northchannel_station_list_amp, northchannel_station_list_phase,...
    output_pred_disp)

%{
locations = cell(size(stationlist));
locations{1} = struct('station','NPB', 'latlon', [19.41200 -155.28100 1115.0],'misfit_weight',1);
locations{2} = struct('station','SRM', 'latlon', [19.39531 -155.27400 1128.0],'misfit_weight',1);
locations{3} = struct('station','NPT', 'latlon', [19.41200 -155.28100 1115.0],'misfit_weight',1);
locations{4} = struct('station','OBL', 'latlon', [19.41750 -155.28400 1107.0],'misfit_weight',1);
locations{5} = struct('station','WRM', 'latlon', [19.40650 -155.30000 1163.0],'misfit_weight',1);
locations{6} = struct('station','SDH', 'latlon', [19.38950 -155.29400 1133.0],'misfit_weight',1);
locations{7} = struct('station','UWE', 'latlon', [19.42111 -155.29300 1240.0],'misfit_weight',1);
locations{8} = struct('station','UWB', 'latlon', [19.42469 -155.27800 1091.0],'misfit_weight',1);
locations{9} = struct('station','SBL', 'latlon', [19.42700 -155.26800 1084.0],'misfit_weight',1);
locations{10} = struct('station','HAT', 'latlon', [19.42300 -155.26100 1082.0],'misfit_weight',0);
locations{11} = struct('station','BYL', 'latlon', [19.41200 -155.26000 1059.0],'misfit_weight',0);
locations{12} = struct('station','KKO', 'latlon', [19.39781 -155.26600 1146.0], 'misfit_weight',1);
locations{13} = struct('station','RIMD', 'latlon', [19.39531 -155.27400 1128.0],'misfit_weight',1);
locations{14} = struct('station','PUHI', 'latlon', [19.385510 -155.251310 1079],'misfit_weight',1);
%convert lat-lon locations to UTM E-N coordinates
zone = utmzone(locations{1}.latlon(1),locations{1}.latlon(2));
utmstruct = defaultm('utm');
utmstruct.zone = zone;
utmstruct.geoid = wgs84Ellipsoid;
utmstruct = defaultm(utmstruct);
for i = 1:length(locations)
    locations{i}.loc = [0 0 0];
   [locations{i}.loc(1), locations{i}.loc(2), locations{i}.loc(3)] = mfwdtran(utmstruct,locations{i}.latlon(1),locations{i}.latlon(2),locations{i}.latlon(3)); 
end
%}

%station data
VariableNames = TQ_tab.Properties.VariableNames;
stationcell = {};
for ii = 1:length(VariableNames)
    if length(VariableNames{ii}) > 7
        if strcmp(VariableNames{ii}(end-3:end), '_HHZ')
            if strcmp(VariableNames{ii}(1:4), 'amp_')
                stationcell{end+1} = VariableNames{ii}(8:10);
            end
        end
    end
end
stations = struct('station',stationcell);

%station locations
for ii = 1:length(stations)
    if strcmp(stations(ii).station,'NPB')
        stations(ii).latlon = [19.41200 -155.28100 1115.0];
    end
    if strcmp(stations(ii).station,'SRM')
        stations(ii).latlon = [19.39531 -155.27400 1128.0];
    end
    if strcmp(stations(ii).station,'NPT')
        stations(ii).latlon = [19.41200 -155.28100 1115.0];
    end
    if strcmp(stations(ii).station,'OBL')
        stations(ii).latlon = [19.41750 -155.28400 1107.0];
    end
    if strcmp(stations(ii).station,'WRM')
        stations(ii).latlon = [19.40650 -155.30000 1163.0];
    end
    if strcmp(stations(ii).station,'SDH')
        stations(ii).latlon = [19.38950 -155.29400 1133.0];
    end
    if strcmp(stations(ii).station,'UWE')
        stations(ii).latlon = [19.42111 -155.29300 1240.0];
    end
    if strcmp(stations(ii).station,'UWB')
        stations(ii).latlon = [19.42469 -155.27800 1091.0];
    end
    if strcmp(stations(ii).station,'SBL')
        stations(ii).latlon = [19.42700 -155.26800 1084.0];
    end
    if strcmp(stations(ii).station,'HAT')
        stations(ii).latlon = [19.42300 -155.26100 1082.0];
    end
    if strcmp(stations(ii).station,'BYL')
        stations(ii).latlon = [19.41200 -155.26000 1059.0];
    end
    if strcmp(stations(ii).station,'KKO')
        stations(ii).latlon = [19.39781 -155.26600 1146.0];
    end
    if strcmp(stations(ii).station,'RIM')
        stations(ii).latlon = [19.39531 -155.27400 1128.0];
        stations(ii).station = 'RIMD';
    end
    if strcmp(stations(ii).station,'PUH')
        stations(ii).latlon = [19.385510 -155.251310 1079];
        stations(ii).station = 'PUHI';
    end
end
%convert lat-lon locations to UTM E-N coordinates
zone = utmzone(stations(ii).latlon(1),stations(ii).latlon(2));
utmstruct = defaultm('utm');
utmstruct.zone = zone;
utmstruct.geoid = wgs84Ellipsoid;
utmstruct = defaultm(utmstruct);
for ii = 1:length(stations)
   [stations(ii).locE, stations(ii).locN, stations(ii).locZ] = ...
       mfwdtran(utmstruct,stations(ii).latlon(1),stations(ii).latlon(2),stations(ii).latlon(3)); 
end

obs = zeros(3,length(stations)); %observation (station) locations
obstilt = zeros(3,4*length(stations));%observation locations for tilt
for i = 1:length(stations)
    obs(1,i) = stations(i).locE;
    obs(2,i) = stations(i).locN;
    obs(3,i) = 0;
    obstilt(:,4*(i-1)+1) = obs(:,i) + [-1;0;0];
    obstilt(:,4*(i-1)+2) = obs(:,i) + [1;0;0];
    obstilt(:,4*(i-1)+3) = obs(:,i) + [0;-1;0];
    obstilt(:,4*(i-1)+4) = obs(:,i) + [0;1;0];
end 

%lava lake location (from Chouet work)
epilatlon = [19.405402 -155.281041 100.0];
[epiloc_east, epiloc_north, ~] = mfwdtran(utmstruct,epilatlon(1),epilatlon(2),epilatlon(3));
mu = 10*10^9; %shear modulus
nu = 0.25; %poisson's ratio
%mogi location (from Kyle Anderson Geodetic inversions for 2018 eruptions)
kyle_east = 63.7 + 260740;
eastoffset = kyle_east-epiloc_east;
kyle_north = -41.95 + 2147480;
northoffset = kyle_north-epiloc_north;
centroideast = eastoffset+epiloc_east;
centroidnorth = northoffset+epiloc_north;


depth_vec = linspace(500,2500,201);
G_cell = cell(size(depth_vec));
Gtilt_cell = cell(size(depth_vec));
for ii = 1:length(depth_vec)
    
    G_cell{ii} = zeros(3*length(stations),6); %"Green's function" displacement matrix
    Gtilt_cell{ii} = zeros(3*length(stations),6); %"Green's function" tilt anlge matrix

    %calculate Greens functions for 11
    mdl = [depth_vec(ii) centroideast centroidnorth 10^6 0 0 0 0 0]';
    U = quasistatic_MT_greensfunc(mdl,obs,mu,nu);
    Utilt = quasistatic_MT_greensfunc(mdl,obstilt,mu,nu);
    G_cell{ii}(:,1) = [U(1,:).';U(2,:).';U(3,:).'];
    Gtilt_cell{ii}(:,1) = [atan((Utilt(3,2:4:end)-Utilt(3,1:4:end))/2).'; atan((Utilt(3,4:4:end)-Utilt(3,3:4:end))/2).';zeros(length(stations),1)];

    %calculate greens functions for 22
    mdl = [depth_vec(ii) centroideast centroidnorth 0 10^6 0 0 0 0]';
    U = quasistatic_MT_greensfunc(mdl,obs,mu,nu);
    Utilt = quasistatic_MT_greensfunc(mdl,obstilt,mu,nu);
    G_cell{ii}(:,2) = [U(1,:).';U(2,:).';U(3,:).'];
    Gtilt_cell{ii}(:,2) = [atan((Utilt(3,2:4:end)-Utilt(3,1:4:end))/2).'; atan((Utilt(3,4:4:end)-Utilt(3,3:4:end))/2).';zeros(length(stations),1)];

    %calculate greens functions for 33
    mdl = [depth_vec(ii) centroideast centroidnorth 0 0 10^6 0 0 0]';
    U = quasistatic_MT_greensfunc(mdl,obs,mu,nu);
    Utilt = quasistatic_MT_greensfunc(mdl,obstilt,mu,nu);
    G_cell{ii}(:,3) = [U(1,:).';U(2,:).';U(3,:).'];
    Gtilt_cell{ii}(:,3) = [atan((Utilt(3,2:4:end)-Utilt(3,1:4:end))/2).'; atan((Utilt(3,4:4:end)-Utilt(3,3:4:end))/2).';zeros(length(stations),1)];

    %calculate Greens functions for 12
    mdl = [depth_vec(ii) centroideast centroidnorth 0 0 0 10^6 0 0]';
    U = quasistatic_MT_greensfunc(mdl,obs,mu,nu);
    Utilt = quasistatic_MT_greensfunc(mdl,obstilt,mu,nu);
    G_cell{ii}(:,4) = [U(1,:).';U(2,:).';U(3,:).'];
    Gtilt_cell{ii}(:,4) = [atan((Utilt(3,2:4:end)-Utilt(3,1:4:end))/2).'; atan((Utilt(3,4:4:end)-Utilt(3,3:4:end))/2).';zeros(length(stations),1)];

    %calculate Greens functions for 13
    mdl = [depth_vec(ii) centroideast centroidnorth 0 0 0 0 10^6 0]';
    U = quasistatic_MT_greensfunc(mdl,obs,mu,nu);
    Utilt = quasistatic_MT_greensfunc(mdl,obstilt,mu,nu);
    G_cell{ii}(:,5) = [U(1,:).';U(2,:).';U(3,:).'];
    Gtilt_cell{ii}(:,5) = [atan((Utilt(3,2:4:end)-Utilt(3,1:4:end))/2).'; atan((Utilt(3,4:4:end)-Utilt(3,3:4:end))/2).';zeros(length(stations),1)];

    %calculate Greens functions for 23
    mdl = [depth_vec(ii) centroideast centroidnorth 0 0 0 0 0 10^6]';
    U = quasistatic_MT_greensfunc(mdl,obs,mu,nu);
    Utilt = quasistatic_MT_greensfunc(mdl,obstilt,mu,nu);
    G_cell{ii}(:,6) = [U(1,:).';U(2,:).';U(3,:).'];
    Gtilt_cell{ii}(:,6) = [atan((Utilt(3,2:4:end)-Utilt(3,1:4:end))/2).'; atan((Utilt(3,4:4:end)-Utilt(3,3:4:end))/2).';zeros(length(stations),1)];
    
end

g = 9.807; %leave positive

%loop over events
% total_stationlist = [eastchannel_station_str(4:6)'; ...
%     northchannel_station_str(4:6)';...
%     vertchannel_station_str(4:6)'];
% total_channellist = [eastchannel_station_str(end-3:end)'; ...
%     northchannel_station_str(end-3:end)';...
%     vertchannel_station_str(end-3:end)'];
if output_pred_disp
   for ii = 1:length(vertchannel_station_list_amp)
       TQ_tab.(strcat(vertchannel_station_list_amp{ii},'_mtdepth')) = NaN(height(TQ_tab),1);
   end
   for ii = 1:length(vertchannel_station_list_phase)
       TQ_tab.(strcat(vertchannel_station_list_phase{ii},'_mtdepth')) = NaN(height(TQ_tab),1);
   end
   for ii = 1:length(eastchannel_station_list_amp)
       TQ_tab.(strcat(eastchannel_station_list_amp{ii},'_mtdepth')) = NaN(height(TQ_tab),1);
   end
   for ii = 1:length(eastchannel_station_list_phase)
       TQ_tab.(strcat(eastchannel_station_list_phase{ii},'_mtdepth')) = NaN(height(TQ_tab),1);
   end
   for ii = 1:length(northchannel_station_list_amp)
       TQ_tab.(strcat(northchannel_station_list_amp{ii},'_mtdepth')) = NaN(height(TQ_tab),1);
   end
   for ii = 1:length(northchannel_station_list_phase)
       TQ_tab.(strcat(northchannel_station_list_phase{ii},'_mtdepth')) = NaN(height(TQ_tab),1);
   end
end
mt_comp_vec = NaN(length(depth_vec),6);
misfit_vec = NaN(size(depth_vec));
for ii = 1:height(TQ_tab)
    %vector of displacements
    FU_real = [TQ_tab{ii,eastchannel_station_list_amp}' ...
        .*exp(1i*TQ_tab{ii,eastchannel_station_list_phase}');...
        TQ_tab{ii,northchannel_station_list_amp}' ...
        .*exp(1i*TQ_tab{ii,northchannel_station_list_phase}');...
        TQ_tab{ii,vertchannel_station_list_amp}' ...
        .*exp(1i*TQ_tab{ii,vertchannel_station_list_phase}')];
    FU_real = FU_real./(2*pi*1i*(1/TQ_tab.T(ii))); %convert from vel to disp
    %remove NaNs
    nanind = isnan(abs(FU_real));
    if sum(~nanind) > 6 %enough channels to invert for
        FU_real = FU_real(~nanind);

        for jj = 1:length(depth_vec)

            G_local = G_cell{jj}(~nanind,:);
            Gtilt_local = Gtilt_cell{jj}(~nanind,:);
        %     stationlist_local = total_stationlist(nanind);
        %     channellist_local = total_stationlist(nanind);

            %invert for deltaV
            mt_comp_vec(jj,:) = ((G_local + g*Gtilt_local/(2*pi*1i*(1/TQ_tab.T(ii)))^2)\FU_real)';
            %use to get predicted disp
            FU_pred = (G_local + g*Gtilt_local/(2*pi*1i*(1/TQ_tab.T(ii)))^2)*mt_comp_vec(jj,:)';
            %calculate misfit as total error/total energy
            misfit_vec(jj) = sum(abs(FU_real-FU_pred))/sum(abs(FU_real));
        end
        [~,min_misfit_in] = min(misfit_vec);
        misfit = misfit_vec(min_misfit_in);
        mt_comp = mt_comp_vec(min_misfit_in,:);
        %MT structure for decomp calc
        t = linspace(0,2*pi,40)';
        MT = table();
        MT.MT11 = abs(mt_comp(1))*cos(angle(mt_comp(1)) + t);
        MT.MT22 = abs(mt_comp(2))*cos(angle(mt_comp(2)) + t);
        MT.MT33 = abs(mt_comp(3))*cos(angle(mt_comp(3)) + t);
        MT.MT12 = abs(mt_comp(4))*cos(angle(mt_comp(4)) + t);
        MT.MT13 = abs(mt_comp(5))*cos(angle(mt_comp(5)) + t);
        MT.MT23 = abs(mt_comp(6))*cos(angle(mt_comp(6)) + t);
        MT = MT_decomp(MT);
        TQ_tab.mtdepth_misfit(ii) = misfit;
        TQ_tab.mtdepth_depth(ii) = depth_vec(min_misfit_in);
        TQ_tab.mtdepth_MT11(ii) = mt_comp(1);
        TQ_tab.mtdepth_MT22(ii) = mt_comp(2);
        TQ_tab.mtdepth_MT33(ii) = mt_comp(3);
        TQ_tab.mtdepth_MT12(ii) = mt_comp(4);
        TQ_tab.mtdepth_MT13(ii) = mt_comp(5);
        TQ_tab.mtdepth_MT21(ii) = mt_comp(6);
        TQ_tab.mtdepth_iso(ii) = mean(abs(MT.iso));
        TQ_tab.mtdepth_clvd(ii) = mean(abs(MT.clvd));
        TQ_tab.mtdepth_dc(ii) = mean(abs(MT.dc));
        TQ_tab.mtdepth_scalarmoment(ii) = mean(abs(MT.scalarmoment));
    %     misfit = sum(abs(sqrt((abs(FU_pred).*cos(angle(FU_pred)) - abs(FU_real).*cos(angle(FU_real))).^2 ...
    %         + (-abs(FU_pred).*sin(angle(FU_pred)) + abs(FU_real).*sin(angle(FU_real))).^2))) ...
    %         /sum(abs(FU_real));
    %     TQ_tab.mogi_misfit(ii) = misfit;

        if output_pred_disp
            %add predicted vel to table
            FU_predvel = FU_pred.*(2*pi*1i*(1/TQ_tab.T(ii))); %convert from vel to disp
            FU_pred_nan = NaN(size(nanind));
            nancount = 0;
            for jj = 1:length(nanind)
                if ~nanind(jj)
                    nancount = nancount+1;
                    FU_pred_nan(jj) = FU_predvel(nancount);
                end
            end
            eaststart = 0;
            for jj = 1:length(eastchannel_station_list_amp)
                TQ_tab.(strcat(eastchannel_station_list_amp{jj},'_mtdepth'))(ii) = ...
                    abs(FU_pred_nan(eaststart + jj));
                TQ_tab.(strcat(eastchannel_station_list_phase{jj},'_mtdepth'))(ii) = ...
                    angle(FU_pred_nan(eaststart + jj));
            end
            northstart = length(eastchannel_station_list_amp);
            for jj = 1:length(northchannel_station_list_amp)
                TQ_tab.(strcat(northchannel_station_list_amp{jj},'_mtdepth'))(ii) = ...
                    abs(FU_pred_nan(northstart + jj));
                TQ_tab.(strcat(northchannel_station_list_phase{jj},'_mtdepth'))(ii) = ...
                    angle(FU_pred_nan(northstart + jj));
            end
            vertstart = length(eastchannel_station_list_amp) + length(northchannel_station_list_amp);
            for jj = 1:length(vertchannel_station_list_amp)
                TQ_tab.(strcat(vertchannel_station_list_amp{jj},'_mtdepth'))(ii) = ...
                    abs(FU_pred_nan(vertstart + jj));
                TQ_tab.(strcat(vertchannel_station_list_phase{jj},'_mtdepth'))(ii) = ...
                    angle(FU_pred_nan(vertstart + jj));
            end
        end

        %calculate angles
        count = 0;
        test_time = linspace(0,TQ_tab.T(ii),60);
        keepind = 1:length(nanind);
        keepind = keepind(~nanind);
        for jj = 1:length(stationlist)
            aeastind = jj;
            anorthind = length(stationlist) + jj;
            avertind = 2*length(stationlist) + jj;
            if sum(nanind([aeastind, anorthind, avertind])) == 0 %all channels valid
                count = count+1;
                eastind = find(aeastind == keepind);
                northind = find(anorthind == keepind);
                vertind = find(avertind == keepind);
                %approximate mean angle error numerically (could do algebraically later)
                % by adding up over one period
                east = abs(FU_real(eastind))*cos(1/TQ_tab.T(ii)*(2*pi)*test_time + angle(FU_real(eastind)));
                eastpred = abs(FU_pred(eastind))*cos(1/TQ_tab.T(ii)*(2*pi)*test_time + angle(FU_pred(eastind)));
                north = abs(FU_real(northind))*cos(1/TQ_tab.T(ii)*(2*pi)*test_time + angle(FU_real(northind)));
                northpred = abs(FU_pred(northind))*cos(1/TQ_tab.T(ii)*(2*pi)*test_time + angle(FU_pred(northind)));
                vert = abs(FU_real(vertind))*cos(1/TQ_tab.T(ii)*(2*pi)*test_time + angle(FU_real(vertind)));
                vertpred = abs(FU_pred(vertind))*cos(1/TQ_tab.T(ii)*(2*pi)*test_time + angle(FU_pred(vertind)));
                %calculate angles
                angles = (east.*eastpred + north.*northpred + vert.*vertpred);
                mag = (east.^2 + north.^2 + vert.^2).^0.5;
                magpred = (eastpred.^2 + northpred.^2 + vertpred.^2).^0.5;
                angles = acos(angles./(mag.*magpred));
                if count == 1
                    meanangles = mean(abs(wrapToPi(angles)));
                else
                    meanangles = meanangles + mean(abs(wrapToPi(angles)));
                end
            end
        end
        if count > 0
            meanangles = meanangles/count;
            TQ_tab.mtdepth_mean_angle_misfit(ii) = meanangles*360/(2*pi);
        end
        
    end %end if enough channels
end %end loop over events

end

